#ifndef FIELD
#define FIELD

#include <vector>
#include "glm/glm.hpp"

template <typename T>
class Field {
public:
	enum InterpolationMode { 
		BILIN = 0,
		BICUBIC,
		NUM_INTERP_MODES
	};

	void SetInterpMode(InterpolationMode inMode ){
		mInterpMode =inMode;
	}

	Field() :
	  mDimX(0), mDimY(0)
	  {}

	Field(int w, int h, InterpolationMode im = BILIN) :
		mDimX(w), mDimY(h), mInterpMode(im)
	  {
		  // Allocate vector
		  mData.resize(w*h);
	  }
	
	Field(const Field & copy) 
	{
		(*this) = copy;
	}

	Field(std::vector<T> & inData, int w, int h, InterpolationMode im = BILIN) :
		mDimX(w), mDimY(h), mInterpMode(im)
	  {
		  // Allocate vector
		  mData.resize(w*h);
		  std::copy(inData.begin(), inData.end(), mData.begin());
	  }

	Field & operator=(const Field & rhs){
		if (&rhs != this){
			mData = rhs.mData;
			mDimX = rhs.mDimX;
			mDimY = rhs.mDimY;
			mInterpMode = rhs.mInterpMode;
		}

		return *this;
	}

	T& operator[](unsigned int i){
		return mData[i];
	}

	T& operator()(int i, int j) 
	{
		i = std::min(std::max(i, 0), mDimX-1);
		j = std::min(std::max(j, 0), mDimY-1);
		return mData[ndx(i, j)];
	}

	const T& operator()(int i, int j) const{
		i = std::min(std::max(i, 0), mDimX-1);
		j = std::min(std::max(j, 0), mDimY-1);
		return mData[ndx(i,j)];
	}

	InterpolationMode GetInterpMode() const {
		return mInterpMode;
	}
	
	// Sample via the interpolation mode set.
	void Sample ( float u, float v, T & out ) const {
		switch(GetInterpMode()){
		case BILIN:
			SampleBilinInterp(u, v, out);
			break;
		case BICUBIC:
			SampleBicubicInterp(u, v, out);
			break;
		default:
			break;
		}
	}

	std::vector<T> * GetRawData() {
		return &mData;
	}

	int GetCellU( float u ) const {
		// Convert u into the cell along the x direction.
		return (int)(u*mDimX);
	}

	int GetCellV( float v ) const {
		// Copnvert v into the cell along the y direction.
		return (int)(v*mDimY); 
	}

	unsigned int GetDimX() const {
		return mDimX;
	}

	unsigned int GetDimY() const {
		return mDimY;
	}

private:
	unsigned int ndx( unsigned int i, unsigned int j ) const{
		return j * mDimX + i;
	}

	void SampleBicubicInterp( float u, float v, T & out) const {
		// First get the 4 cell corners we need.
		int cell[2] = {
			GetCellU(u), GetCellV(v) 
		};

		// Convert to be local to the cell its in.
		u = (u*mDimX) - (float)cell[0];
		v = (v*mDimY) - (float)cell[1];

		// Get the velocities at the cell's corners.
		int cOffset[2];

		cOffset[0] = (u <= 0.5f) ? -1 : 1;
		cOffset[1] = (v <= 0.5f) ? -1 : 1;

		//left, right, bottom, top
		glm::vec2 sampleNdx = glm::vec2( (cOffset[0] < 0) ? cell[0]-1 : cell[0], (cOffset[1] < 0 ? cell[1]-1 : cell[1]));

		// Translate over
		float su = u > 0.5 ? u- 0.5f : u+0.5f,
			sv = v > 0.5 ? v - 0.5f : v+0.5f;

		const Field<T> & ref = *this;

		BicubicInterp(
			sampleNdx[0],
			sampleNdx[1],
			ref,
			su, 
			sv,
			out );
	}

	void SampleBilinInterp ( float u, float v, T & out ) const {
		// First get the 4 cell corners we need.
		int cell[2] = {
			GetCellU(u), GetCellV(v) 
		};

		// Convert to be local to the cell its in.
		u = (u*mDimX) - (float)cell[0];
		v = (v*mDimY) - (float)cell[1];

		// Get the velocities at the cell's corners.
		int cOffset[2];

		cOffset[0] = (u <= 0.5f) ? -1 : 1;
		cOffset[1] = (v <= 0.5f) ? -1 : 1;

		//left, right, bottom, top
		int sampleNdx[4][2] = 
		{
			{(cOffset[0] < 0) ? cell[0]-1 : cell[0], (cOffset[1] < 0 ? cell[1]-1 : cell[1]) },
			{(cOffset[0] < 0) ? cell[0] : cell[0] + 1, (cOffset[1] < 0 ? cell[1]-1 : cell[1])},
			{(cOffset[0] < 0) ? cell[0] : cell[0] + 1, (cOffset[1] < 0 ? cell[1] : cell[1]+1) },
			{(cOffset[0] < 0) ? cell[0]-1 : cell[0],(cOffset[1] < 0 ? cell[1] : cell[1]+1)}
		};

		// Translate over
		float su = u > 0.5 ? u- 0.5f : u+0.5f,
			sv = v > 0.5 ? v - 0.5f : v+0.5f;

		const Field & ref = *this;

		BilinInterp(
			ref(sampleNdx[0][0], sampleNdx[0][1]),
			ref(sampleNdx[1][0], sampleNdx[1][1]),
			ref(sampleNdx[2][0], sampleNdx[2][1]),
			ref(sampleNdx[3][0], sampleNdx[3][1]),
			su, 
			sv,
			out );
	}

	void BilinInterp( const T & lowerL,
		const T & lowerR,
		const T & upperR,
		const T & upperL, 
		float u,
		float v,
		T & out ) const {
		
			out = (1-v)*((1-u)*lowerL + u*(lowerR)) 
				+
				(v)*((1-u)*upperL + u*(upperR));
	}

	void BicubicInterp( 
		int i, // lower left cell of interpolation region
		int j,
		const Field<T> & inField,
		float u,
		float v,
		T & out ) const {

			T cuby[4];
			
			CubicInterp(inField(i-1, j-1), inField(i, j-1), inField(i+1, j-1), inField(i+2, j-1), u, cuby[0]);
			CubicInterp(inField(i-1, j), inField(i, j), inField(i+1, j), inField(i+2, j), u, cuby[1]);
			CubicInterp(inField(i-1, j+1), inField(i, j+1), inField(i+1, j+1), inField(i+2, j+1), u, cuby[2]);
			CubicInterp(inField(i-1, j+2), inField(i, j+2), inField(i+1, j+2), inField(i+2, j+2), u, cuby[3]);

			CubicInterp(cuby[0], cuby[1], cuby[2], cuby[3], v, out);
	}

	void CubicInterp( const glm::vec2 & f0, 
						const glm::vec2 & f1,
						const glm::vec2 & f2,
						const glm::vec2 & f3,
						float t,
						glm::vec2 & out ) const{
		float h = 1.0f/(mDimX);
		
		glm::vec2 df1 = (f2-f0)*0.5f,
		  df2 = (f3-f1)*0.5f;

		if ( glm::dot(df1,(f2-f1)) <0){ df1 = glm::vec2(0,0); }
		if ( glm::dot(df2,(f2-f1)) <0){ df2 = glm::vec2(0,0); }

		out = 
			(2.0f*f1 - 2.0f*f2 + df1 + df2) * (t*t*t) +
			(-3.0f*f1 + 3.0f*f2 - 2.0f*df1 - df2) * (t*t) +
			df1*t +
			f1;

		//out = 
		//	(-0.5f*f0 + (3.0f/2.0f)*f1 - (3.0f/2.0f)*f2 + 0.5f*f3)*(t*t*t) +
		//	(f0 - (5.0f/2.0f)*f1 + 2.0f*f2 - 0.5f*f3)*(t*t) +
		//	((-0.5f)*f0 + 0.5f*f2)*t +
		//	f1;
	}

	void CubicInterp( const float & f0, 
						const float & f1,
						const float & f2,
						const float & f3,
						float t,
						float & out ) const{
		float h = 1.0f/(mDimX);
		
		float df1 = (f2-f0)*0.5f,
		  df2 = (f3-f1)*0.5f;

		if ( df1*(f2-f1) <0){ df1 = 0; }
		if ( df2*(f2-f1) <0){df2=0;}

		out = 
			(2.0f*f1 - 2.0f*f2 + df1 + df2) * (t*t*t) +
			(-3.0f*f1 + 3.0f*f2 - 2.0f*df1 - df2) * (t*t) +
			df1*t +
			f1;

		//out = 
		//	(-0.5f*f0 + (3.0f/2.0f)*f1 - (3.0f/2.0f)*f2 + 0.5f*f3)*(t*t*t) +
		//	(f0 - (5.0f/2.0f)*f1 + 2.0f*f2 - 0.5f*f3)*(t*t) +
		//	((-0.5f)*f0 + 0.5f*f2)*t +
		//	f1;
	}


	InterpolationMode mInterpMode;

	std::vector<T> mData;
	int mDimX, mDimY;
};


#endif